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Abstract 



Two different models exhibiting self-organized criticality are analyzed by 
means of the dynamic renormalization group. Although the two models differ 
by their behavior under a parity transformation of the order parameter, it 
is shown that they both belong to the same universality class, in agreement 
with computer simulations. The asymptotic values of the critical exponents 
are estimated up to one loop order from a systematic expansion of a nonlinear 
equation in the number of coupling constants. 
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Recently much attention has been paid to the phenomenon known as self-organized crit- 
icality. Bak, Tang, and Wiesenfeld [|l|] studied a cellular automaton model as a paradigm 
for the explanation of two widely occurring phenomena in nature: 1/f noise and fractal 
structures. Both have in common a lack of characteristic scales. This scale invariance sug- 
gests that the system is critical in analogy with classical equilibrium critical phenomena, but 
in self-organized criticality one deals with dynamical nonequilibrium statistical properties. 
On the other hand the system evolves naturally to the critical state without any tuning of 
external parameters. 

Several cellular automata models exhibiting self-organized criticality have been reported 
in the literature. In the original sandpile model of Bak et al. ]I] the system is perturbed 
externally by a random addition of sand grains. Once the slope between two contiguous 
cells have reached a threshold condition, sand is transferred to its neighbors in a fixed 
amount. Taking this model as a reference, different dynamical rules have been investigated 
leading to a wide variety of universality classes. Some authors have attempted to connect 
these models to stochastic differential equations @,[|. These continuous descriptions are 
built according to the symmetry rules obeyed by the discrete models. Nevertheless none of 
them explicitly includes the threshold condition. On the other hand anomalous diffusion 
equations with singularities in the diffusion coefficient have been considered in order to study 
the deterministic dynamics of the avalanches generated in the critical state 

In this paper we study two nonlinear stochastic differential equations derived from the 
discrete dynamical rules of two models which have different symmetry properties from which 
one expects different critical behavior, but in the simulations they give the same critical 
exponents. We will show analytically, by means of the dynamic renormalization group 
(DRG) [§,[7| , that both models belong to the same universality class. The threshold condition 
is kept but the step function is regularized in order to allow a power series expansion. It 
is in the limit where the threshold is recovered that the coupling constants that distinguish 
both models are shown to be irrelevant. Once this equivalence is established we deal with 
one of the models, showing that an infinite number of coupling constants is relevant below 
the upper critical dimension d c = 4; by expanding in the number of nonlinearities to first 
order in e = 4 — d, we get an accurate estimate of the dynamical exponent close to the value 
obtained in the numerical simulations. 

First, we describe briefly the dynamics of the two models showing self-organized critical- 
ity. The first model, to be referred to as model A, was originally proposed by Zhang || and 
it consists of a lattice in which any site can store some energy E continuously distributed 
between and E c . This variable which we call energy can have different physical interpreta- 
tions ||. The system is perturbed by adding a random amount of energy SE at a randomly 
chosen site. Once a site reaches a value of the energy greater than some critical value E c this 
site becomes active and transfers the energy to its nearest neighbors. At this point the input 
of energy from the outside is turned off and the energy transferred to the neighboring sites 
can make them active giving rise to an activation cluster or avalanche that ends when all 
the sites have reached a value of the energy smaller than E c . It is only when the avalanche 
has stopped that energy is added again, otherwise the system remains quiescent. Model B is 
closer to the original sandpile model of Bak et al. M and differs from A in that the amount 
of energy a critical site transfers to its neighbors is a fixed amount E c instead of its whole 
energy E. 
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The microscopic rules for these models can be written in the form of a set of algebraic 
equations and coarse-grained to get effective equations in terms of a rescaled energy E—E c ^ 
E as follows 

= «V 2 [9(E(v,t)) (E(r,t) + E c )\ +r,(r,t), (la) 

and 



dt 



aV 2 [9(E(r,t))E c ]+ V (r,t) (lb) 



for models A and B, respectively. In these equations a plays the role of a diffusion constant. 
The threshold condition enters through the Heaviside step function 9. The external noise 
has a zero mean value and a correlation function given by 

< ?7(r, t)r](r', t') >= 2T5 d (r - r') (2) 

as was extensively discussed in [fETf] . 

In order to perform a perturbative expansion of Eq. (1) we choose one of the represen- 
tations of the step function |J 



9{x) = lim — 

,3—>oo 7T 



7T 

— + arctan px 



(3) 



and make a series expansion of the arctan in powers of its argument. 
We can then formally write for both models 



dE( > r > *) = DV 2 £(r, t) + f; A n V 2 £"(r, t) + ,(r, t) (4) 



where the effective diffusion constant (D) and the coupling constants (A n ) take different 
values depending on the symmetry of the model 

D* = a(^ + l);D° = ^ (5) 



A {-\YaE c ^ B _ A B {-lYaE c ^ 

A 2n - n{2n _ 1) » A 2n-U, A 2ri _i — A 2ri _i — 7r(2n _ 1) W 

The difference between both models is that B has a definite parity under E — > — 
because for this model the even coupling constants vanish. Moreover even and odd coupling 
constants depend on j3 in a different way. This fact will be shown to be crucial as far as 
universality classes are concerned since for the original models j3 — > oo. In this limit we can 
see that D A = D B = D. 

The scaling behavior of (1) is obtained from the correlation function 

< (E(r , t ) - £(r + r, t + t)f >s~ r^f(t/r z ) (7) 
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where \ and z are the roughening and dynamical exponents, respectively. The relevance of 
the different coupling constants can be checked by naive dimensional analysis: a change of 
scale r — > e l r is followed by t — > e zl t and E — > e xl E in order for (0) be satisfied. Under 
this scaling transformation z and x are chosen to keep the linear equation invariant. With 
these values one can check that all dimensionless coupling constants are of the same order 
in e = 4 — d and hence all nonlinear terms are equally relevant for d < 4. 

The dynamic renormalization group procedure consists of two steps: an elimination of 
the fast short wave-length modes followed by a rescaling of the remaining modes The 
recursion relations after an infinitesimal transformation for the coefficients are, to one loop 
order 
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dD 

~dT 



D 



2 + 3- 



rA d-4 A 



D 3 



TA 

- -4 



(8a) 



dT 



2 X ~d] 



(8b) 



dl 



n+2 



X n [(n-l) X + z- 2)] + £ 7j rA d - 4 Z^ J[ \\ 



i=2 



(8c) 



where jj are numerical factors and by are natural numbers accounting for the number of 
times the bare vertex Aj enters the renormalization of the vertex A„ in the j-th term of 
the sum. By simple dimensional analysis we get aj = — 1 — J^bij and ~ — n + 1 



obtaining that the terms within the sum of 
over even values of i. Thus we obtain 



behave as (3 n So 6 ^ where the sum runs 



1 d\ n 
X n dl 




13 n even 
n odd 



(9) 



And from this behavior we can conclude: a) for odd n any term containing even vertices 
will not contribute when compared with terms involving only odd vertices when (3 — ► oo, 
in particular A3 and D = Ai get decoupled from all even coupling constants; b) for even n 
only those terms with one even vertex will contribute [J In this case we can ensure that 
there are no divergences caused by the j3 — > 00 limit in the equations for the even coupling 
constants and thus no even terms are generated under renormalization in the equations for 
the odd coupling constants. The critical exponents are calculated at the fixed points and 
from the previous arguments they will only depend on A3 which in turn is only coupled to 
odd coupling constants. We can conclude finally that the values of the critical exponents 
will be the same irrespective of the values taken by the even nonlinearities showing thus the 



1 The parity of a renormalized vertex is — (— l) e where e is the number of even bare vertices. Thus 
an even renormalized vertex always contains an odd number of even bare vertices. 
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equivalence between the Zhang and BTW models. After having shown that both models 
belong to the same universality class we will focus our analysis on simpler model B. 

Our starting point is Eq. (|J) with X 2n = 0. Since all terms in this expansion are equally 
relevant below the upper critical dimension (d c = 4) a calculation of the critical exponents 
would involve all the contributions. Nevertheless the exponents can be obtained to first 
order in e as a series in the number of coupling constants which are taken into account. 

One starts with the simpler case 



dE{r,t) 



Of 



where D = ^ anc i a, 



= DV 2 E(r } t) + X 3 V A E 6 (r, t) + V (r, t) 

aE c j3 3 



(10) 



3n 



By application of the DRG procedure we get the following differential equation for the 

rA d ~ 4 A 3 

L>3 



effective dimensionless coupling constant Aq 







dl 



e - 27A, 



(11) 



subjected to the initial condition A° = — ^2^2 4 • From ( |TTD we get the following behavior 
in A3 space. For d < 4 one can expect different behavior depending on the initial sign of 
A3. When it is positive it flows towards a stable fixed point A3 = e/27 giving a dynamical 
exponent z = 2 — e/9. A negative A3 flows away giving rise in principle to a strong coupling 
fixed point. The exact solution of Eq. (|TT|) shows that A3, even for A3 < 0, flows back to 
the stable fixed point. But in this case the physical relevance of the stable fixed point is 
a conjecture since along the flow in parameter space the effective coupling constant should 
stay of order e. On the other hand, above the upper critical dimension the system evolves 
towards the trivial fixed point A3 = giving a diffusive behavior with z = 2. 

The next step is to study the contribution of A5 to the critical exponents. The equation 
to study is equivalent to ( PH ) up to fifth order. In this case the differential recursion relations 
for the effective coupling constants 



A, 



rA d - 4 A 3 

D 3 



A, 



rA d - 4 A 5 

£> 2 A.3 



(12) 



are |TT| 



dh 
dl 



A, 



e - 27A 3 + 10A E 



and 



dX 5 
dl 



78A, 



A 2 

lOAs + 81^ 
As 



(13) 



with initial conditions A3 



■TA d 



3a 2 E? 



and A° 



37r 2 rA d - 4 

5a 2 E? 



For d < 4 and A3 > the solution of these equations flows towards a new stable fixed 
point whose projection onto the A3 axis is shifted with respect to the one calculated only 
with the third-order nonlinearity, giving a different dynamical exponent 



159 + 3V1009 
1350 



-e. 



(14) 
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Nevertheless for our model with negative initial values of A3 and A5 the effective coupling 
constants flow away to a strong coupling fixed point under RG transformations. The same 
line of reasoning as for A3 can be applied to the present case and the values of the critical 
exponents can be obtained from the stable fixed point. 

We have performed calculations with A7 and Ag which give further corrections to the 
dynamical exponent |T|. The results are presented in Fig. 1. where one can see that the 



asymptotic behavior approaches the values obtained in M by means of scaling arguments 
and to the results from computer simulations ||10|| . The agreement is much better in d = 3 
than in d — 2 since we are performing a one-loop expansion. Up to the above mentioned 
order we have found that all possible choices of the dimensionless coupling constants have an 
initial value which is independent of f3, legitimizing thus the expansion of the step function 

To summarize, we have studied the critical exponents of two models which show self- 
organized criticality. From the microscopic rules one writes effective long wave-length equa- 
tions involving the threshold condition. After a power series expansion of the threshold 
condition, these equations are suitable for an application of the DRG. This allows us to 
show that both models are in the same universality class and to get an estimate for the 
dynamical exponents to order e by expanding in the number of coupling constants. The 
application of this technique should also be useful for other kinds of problems in which an 
infinite number of relevant nonlinearities arises. 

The author wishes to acknowledge J. Bonet Avalos for very fruitful discussions and A.- 
M.S. Tremblay for a critical reading of the manuscript. This work has been supported by 
CICyT of the Spanish Government, grants #PB89-0233 and #PB92-0863. 
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FIGURES 

FIG. 1. Dynamical exponent as a function of 1/n, where n is the highest odd vertex used to 
compute the dynamical exponent to one-loop order. Circles and squares are the analytical results, 
dashed lines are best linear fits, and the crosses are the results from [8], z = 



8 



